home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 1998 November / IRIX 6.5.2 Base Documentation November 1998.img / usr / share / catman / p_man / cat3 / complib / zlatbs.z / zlatbs
Text File  |  1998-10-30  |  8KB  |  265 lines
  1.  
  2.  
  3.  
  4. ZZZZLLLLAAAATTTTBBBBSSSS((((3333FFFF))))                                                          ZZZZLLLLAAAATTTTBBBBSSSS((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      ZLATBS - solve one of the triangular systems   A * x = s*b, A**T * x =
  10.      s*b, or A**H * x = s*b,
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE,
  14.                         CNORM, INFO )
  15.  
  16.          CHARACTER      DIAG, NORMIN, TRANS, UPLO
  17.  
  18.          INTEGER        INFO, KD, LDAB, N
  19.  
  20.          DOUBLE         PRECISION SCALE
  21.  
  22.          DOUBLE         PRECISION CNORM( * )
  23.  
  24.          COMPLEX*16     AB( LDAB, * ), X( * )
  25.  
  26. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  27.      ZLATBS solves one of the triangular systems
  28.  
  29.      with scaling to prevent overflow, where A is an upper or lower triangular
  30.      band matrix.  Here A' denotes the transpose of A, x and b are n-element
  31.      vectors, and s is a scaling factor, usually less than or equal to 1,
  32.      chosen so that the components of x will be less than the overflow
  33.      threshold.  If the unscaled problem will not cause overflow, the Level 2
  34.      BLAS routine ZTBSV is called.  If the matrix A is singular (A(j,j) = 0
  35.      for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is
  36.      returned.
  37.  
  38.  
  39. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  40.      UPLO    (input) CHARACTER*1
  41.              Specifies whether the matrix A is upper or lower triangular.  =
  42.              'U':  Upper triangular
  43.              = 'L':  Lower triangular
  44.  
  45.      TRANS   (input) CHARACTER*1
  46.              Specifies the operation applied to A.  = 'N':  Solve A * x = s*b
  47.              (No transpose)
  48.              = 'T':  Solve A**T * x = s*b  (Transpose)
  49.              = 'C':  Solve A**H * x = s*b  (Conjugate transpose)
  50.  
  51.      DIAG    (input) CHARACTER*1
  52.              Specifies whether or not the matrix A is unit triangular.  = 'N':
  53.              Non-unit triangular
  54.              = 'U':  Unit triangular
  55.  
  56.      NORMIN  (input) CHARACTER*1
  57.              Specifies whether CNORM has been set or not.  = 'Y':  CNORM
  58.              contains the column norms on entry
  59.              = 'N':  CNORM is not set on entry.  On exit, the norms will be
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. ZZZZLLLLAAAATTTTBBBBSSSS((((3333FFFF))))                                                          ZZZZLLLLAAAATTTTBBBBSSSS((((3333FFFF))))
  71.  
  72.  
  73.  
  74.              computed and stored in CNORM.
  75.  
  76.      N       (input) INTEGER
  77.              The order of the matrix A.  N >= 0.
  78.  
  79.      KD      (input) INTEGER
  80.              The number of subdiagonals or superdiagonals in the triangular
  81.              matrix A.  KD >= 0.
  82.  
  83.      AB      (input) COMPLEX*16 array, dimension (LDAB,N)
  84.              The upper or lower triangular band matrix A, stored in the first
  85.              KD+1 rows of the array. The j-th column of A is stored in the j-
  86.              th column of the array AB as follows:  if UPLO = 'U', AB(kd+1+i-
  87.              j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j)
  88.              = A(i,j) for j<=i<=min(n,j+kd).
  89.  
  90.      LDAB    (input) INTEGER
  91.              The leading dimension of the array AB.  LDAB >= KD+1.
  92.  
  93.      X       (input/output) COMPLEX*16 array, dimension (N)
  94.              On entry, the right hand side b of the triangular system.  On
  95.              exit, X is overwritten by the solution vector x.
  96.  
  97.      SCALE   (output) DOUBLE PRECISION
  98.              The scaling factor s for the triangular system A * x = s*b,  A**T
  99.              * x = s*b,  or  A**H * x = s*b.  If SCALE = 0, the matrix A is
  100.              singular or badly scaled, and the vector x is an exact or
  101.              approximate solution to A*x = 0.
  102.  
  103.      CNORM   (input or output) DOUBLE PRECISION array, dimension (N)
  104.  
  105.              If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains
  106.              the norm of the off-diagonal part of the j-th column of A.  If
  107.              TRANS = 'N', CNORM(j) must be greater than or equal to the
  108.              infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be
  109.              greater than or equal to the 1-norm.
  110.  
  111.              If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns
  112.              the 1-norm of the offdiagonal part of the j-th column of A.
  113.  
  114.      INFO    (output) INTEGER
  115.              = 0:  successful exit
  116.              < 0:  if INFO = -k, the k-th argument had an illegal value
  117.  
  118. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  119.      A rough bound on x is computed; if that is less than overflow, ZTBSV is
  120.      called, otherwise, specific code is used which checks for possible
  121.      overflow or divide-by-zero at every operation.
  122.  
  123.      A columnwise scheme is used for solving A*x = b.  The basic algorithm if
  124.      A is lower triangular is
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. ZZZZLLLLAAAATTTTBBBBSSSS((((3333FFFF))))                                                          ZZZZLLLLAAAATTTTBBBBSSSS((((3333FFFF))))
  137.  
  138.  
  139.  
  140.           x[1:n] := b[1:n]
  141.           for j = 1, ..., n
  142.                x(j) := x(j) / A(j,j)
  143.                x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  144.           end
  145.  
  146.      Define bounds on the components of x after j iterations of the loop:
  147.         M(j) = bound on x[1:j]
  148.         G(j) = bound on x[j+1:n]
  149.      Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  150.  
  151.      Then for iteration j+1 we have
  152.         M(j+1) <= G(j) / | A(j+1,j+1) |
  153.         G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  154.                <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  155.  
  156.      where CNORM(j+1) is greater than or equal to the infinity-norm of column
  157.      j+1 of A, not counting the diagonal.  Hence
  158.  
  159.         G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  160.                      1<=i<=j
  161.      and
  162.  
  163.         |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  164.                                       1<=i< j
  165.  
  166.      Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTBSV if the
  167.      reciprocal of the largest M(j), j=1,..,n, is larger than
  168.      max(underflow, 1/overflow).
  169.  
  170.      The bound on x(j) is also used to determine when a step in the columnwise
  171.      method can be performed without fear of overflow.  If the computed bound
  172.      is greater than a large constant, x is scaled to prevent overflow, but if
  173.      the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a
  174.      non-trivial solution to A*x = 0 is found.
  175.  
  176.      Similarly, a row-wise scheme is used to solve A**T *x = b  or A**H *x =
  177.      b.  The basic algorithm for A upper triangular is
  178.  
  179.           for j = 1, ..., n
  180.                x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  181.           end
  182.  
  183.      We simultaneously compute two bounds
  184.           G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  185.           M(j) = bound on x(i), 1<=i<=j
  186.  
  187.      The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add
  188.      the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.  Then the
  189.      bound on x(j) is
  190.  
  191.           M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. ZZZZLLLLAAAATTTTBBBBSSSS((((3333FFFF))))                                                          ZZZZLLLLAAAATTTTBBBBSSSS((((3333FFFF))))
  203.  
  204.  
  205.  
  206.                <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  207.                          1<=i<=j
  208.  
  209.      and we can safely call ZTBSV if 1/M(n) and 1/G(n) are both greater than
  210.      max(underflow, 1/overflow).
  211.  
  212.  
  213.  
  214.  
  215.  
  216.  
  217.  
  218.  
  219.  
  220.  
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.